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We have developed a method to simulate multiple electron scattering in a vacuum barrier using 
real-space single-electron wavefunctions of the separate surfaces. The tunnelling current is calculated 
to first order in the Dyson series. In zero order, we find a result which differs from the usual Bardeen 
approach by the decay constants of the wavefunctions, entering the description as individual weights 
of tunnelling transitions. To first order we find multiple electron scattering, which can be formulated 
in terms of Bardeen matrices. Here, we also derive a first-principles formulation for the interaction 
energy between the two surfaces. With this method the tunnelling current can in principle be 
computed to any order in the Dyson expansion. 

PACS numbers: 72.10.Bg,72.15.Eb,73.23.Hk 
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From a theoretical point of view a tunnelling electron, 
e.g. in a scanning tunnelling microscopy measurement, 
is part of a system comprising two infinite metal leads 
and an interface, consisting of a vacuum barrier and, op- 
tionally, a molecule or a cluster of atoms with different 
properties than the infinite leads. The system can be 
said to be open - the number of charge carriers is not 
constant - and out of equilibrium - the applied potential 
and charge transport itself introduce polarizations and 
excitations within the system. The theoretical descrip- 
tion of such a system has advanced significantly over the 
last years, to date the most comprehensive description is 
based either on a self-consistent solution of the Lippman- 
Schwinger equation [llor on the non-equilibrium Green's 
function approach Q, H, H, |j, IS - Inelastic effects within 
e.g. a molecule-surface interface can be included by con- 
sidering multiple electron paths from the vacuum into 
the surface substrate [3|- Within the vacuum barrier it- 
self, inelastic effects play an insignificant role. Here, as in 
most experiments in scanning tunneling microscopy, the 
problem can be reduced to the description of the tun- 
nelling current between two leads - the surface 5* and 
the tip T - thought to be in thermal equilibrium. The 
bias potential of the circuit is in this case described by 
a modification of the chemical potentials of surface and 
tip system, symbolized by fis and /zt- This reduces the 
tunnelling problem to the Landauer-Biittiker formulation 
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/ denotes the Fermi distribution functions, G^-^"^^{E) 
and G^{E) are the retarded (advanced) Greens functions 
of the barrier, the Tgi'j-\ are the surface and tip contacts. 
They correspond to the difference of retarded and ad- 
vanced self energy terms of surface and tip; we define 



them by their relation to the spectral function ^5(7^) of 
the surface (tip) Q: 
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Here, the explicit energy dependency has been omit- 
ted for brevity. At present, these equations are evaluated 
within localized basis sets, and in a matrix representa- 
tion. From a theoretical point of view this requires to 
either represent the electronic properties of the two sur- 
faces also in a localized representation 5, 6] , or to trans- 
form the plane wave basis set of most density functional 
methods to a local basis. The use of local basis sets com- 
promises the numerical accuracy in the tunnelling bar- 
rier, since the vacuum tails of the surface wavefunctions 
decay too rapidly: the constant current contours in this 
case are too close to the surface. 

In this Letter we present a formulation of the problem 
which is based on the Green's functions of the two sur- 
faces, given in a real space representation based on the 
electronic eigenstates of the two systems. We show how 
the multiple scattering formalism described in Eq. ^can 
be evaluated in real space, and how it relates to the per- 
turbation expansion of the tunnelling problem. We start 
with an eigenvector expansion of the two surface Green's 
functions, given by: 
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Throughout this paper the wavefunctions ^ and x de- 
note the Kohn-Sham states of surface and tip, respec- 
tively, resulting from a density functional calculation. 
The spectral function As describes the charge density 
matrix, from Eq. [21 we find: 
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The spectral function is related to Ts by Eq. |2| With 
the ansatz for Tg: 
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where C is a constant, we perform the double volume 
integration of Eq. |21 In this case the orthogonality of 
surface states reduces the expression to a compact form: 
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Comparing the result with Eq. |Slwe find for the contacts 
of surface and tip: 

Ts = 2ry^ Vfc(r3)^fe(r4) Tt = 2e^x^(ri)x:(r2) (8) 

k i 

For the construction of the Greens function in the bar- 
rier we use the fact that the charge density is known from 
the separate calculation of surface and tip. In the limit 
of weak coupling, the total charge density of the interface 
is given by: 

n{r,,E) = J2Mrim{ri)S{E~E,) + 

i 

+ ^X,(ri)x;(ri)5(i?-i?,) (9) 
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The surface Green's function fulfills the boundary con- 
ditions at the interface between surface and vacuum. The 
tip Green's function at this position can be considered as 
a weak perturbation, which does not alter the electronic 
properties of the surface. Conversely, the same statement 
holds for the interface between the vacuum and the tip 
with respect to the tip Green's function. This allows to 
construct the zero order approximation for the Green's 
function of the vacuum barrier as the sum of surface and 
tip Green's functions, or: 
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Now all the necessary components for calculating the 
trace in the non equilibrium formalism are given in terms 
of the real space surface and tip wavefunctions. The trace 
involves the following integral: 

Tr [TtG'^TsG'^] = (11) 

= / dri_4rT(ri,r2)G(^)(r2,r3)r5(r3,r4)G'(^)(r4,ri) 

The integrations have to be performed successively. 
We start with the first integral, which involves: 
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Here, we use the physical condition of weak coupling to 
simplify the evaluation. In this limit the overlap integrals 
between states ipm and Xn will be substantially smaller 
than overlaps between Tpm and ipn- Then the integral is 
reduced due to the orthogonality of the wavefunctions to 
the simple form: 
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The same condition leads to the compact result for the 
second integral: 
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The last two integrals involve the overlap between sur- 
face and tip states. The final relation for the trace reads 
then: 
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The calculation of the matrix elements involves an in- 
tegration over infinite space, which cannot directly be 
performed. To convert the volume integrals into surface 
integrals we use the fact that the vacuum states of surface 
and tip are free electron solutions with characteristic de- 
cay constants, complying with the vacuum Schrodinger 
equation. 
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In addition wc make use of the following identities: 

XIV^V-, - V(x:VV;fe) - Vx*V^fe 

^feV^x* = v(V'fcVxn - VxrvV'fe 

After some trivial manipulations, and making use of 
Gauss' theorem, this leads to the result: 
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The surface integral is well known; apart from the uni- 
versal constant h, /2m it describes the tunneling matrix 
element in the perturbation approach jlO. 11] . Denoting 
the surface integral by Mik , we obtain for the zero order 
expansion of the Green's functions the trace: 
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Considering that the states will have discrete energy 
levels E = Ek and setting -q k, e the limit rj.e -^ 0+ 



is easy to evaluate. To check for consistency, we also 
calculated the trace if both contacts are part of the sur- 
face, described by Ts and the Green's function Gs- The 
calculation yields: 
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Considering, that the transition probability in this case 
must be unity, we conclude that the trace has to be scaled 
by a factor of 1/4 to yield the transition probability. The 
final formulation for the tunnelling current thus reads: 
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Note that |A/ifc/(K^ — ^f)!^ is a dimensionless value. 
It describes, like in the Landauer-Biitticker equation [a|, 
the transition probability Tik = tJ}J^ik- Even though the 
formulation looks quite similar to the standard results 
in perturbation theory, it contains some decisive differ- 
ences. While in the perturbation approach the condition 
of resonant tunnelling is explicitly included by a delta 
functional, this is not the case in the scattering formu- 
lation. Instead, the matrix elements are divided by the 
difference between the decay constants of the two wave- 
functions. The decay constants describe the longitudinal 
component of the electron momentum, they enter the de- 
scription, since they account for energy conservation via 
the vacuum Schrodinger equations (see Eq. I16|l . Their 
difference is also related to the speed of charge transfer 
from one side of the junction to the other. Given that the 
denominator differs for every transition, this formulation 
describes, for the first time, how individual contributions 
to the current have to be treated in a proper manner. 

We may compare this new formulation for the zero 
order tunnelling current in the non equilibrium approach 
to the traditional formulation based on Fermi's golden 
rule and the Bardeen method. There we find for the 
current [Ull^: 
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Neglecting initially the decay constants and using 
atomic units [e — m — h — I) the scattering approach 
yields currents for Ei — Ek which are only about one 
tenth the currents obtained in the approach based on 
Fermi's golden rule. However, this is not the full picture. 
The decay constants for metals are in the order of 0.5 to 
0.6 atomic units. Depending on the difference in work 
functions of surface and tip, the current in the scattering 
formulation will be increased by the denominator. The 
increase will be moderate for e.g. tungsten tips and noble 
metal surfaces, where the workfunctions differ by about 



1 to 2 eV. But even in this case it should more than com- 
pensate for the initial difference. The reason seems to be 
that the time of transition is not infinite, as assumed in 
the derivation of Fermi's golden rule, but finite. 

The approach can be extended to higher orders. In 
the first order expansion of the Dyson series the Green's 
function is given by: 
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To calculate the first order Green's function for sys- 
tems out of equilibrium, the equation has to be solved 
self-consistently X,B\. Under tunnelling conditions, how- 
ever, the leads are in thermal equilibrium and the systems 
only weakly coupled. V in this case is the tip potential 
Vt within the vacuum barrier [9j. Self consistency can 
in principle also be achieved by basing the calculation on 
the Kohn-Sham states tp and x of charged surfaces J^j . 

Gfi)(ri,r2) = 

Gfo)(ri,r2)+|dr3Gfo)(ri,r3)yT(r3)Gfo)(r3,r2) (24) 

Setting e ~ rj and with the shortcut /^^ = {E ~ Ei ± 
irj) {E ~ Eki: irj) the first order Green's function depends 
on three additional terms: 
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On the one hand the three additional contributions 
describe excitations of surface electrons. These contri- 
butions can be included in the formulation by a suitable 
adaptation of many-body theory. On the other hand they 
describe the scattering across the barrier due to the vicin- 
ity of surface and tip. These matrix elements, however, 
are well known. Apart from an insignificant contribu- 
tion within the surface system, they are just the Bardeen 
matrix elements in the zero order expansion, or ^ 11] : 
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Limiting the evaluation in the following to scattering 
events across the barrier, we may write for the first order 
Green's function the expansion: 
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The trace now has to be evaluated for the additional 
contributions. The same strategy as used for the zero 



order expansion leads in the limit of r^ ^ 0^ to only one 
surviving multiple scattering term, given by: 
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The current through the tunnelling barrier in the first 
order expansion of the scattering series is therefore: 
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Note that the product in the bottom line is again di- 
mensionless and describes the transition probability in 
the multiple scattering case. The expression contains 
some interesting features. Firstly, the dynamics of elec- 
tron scattering are fully included by the decay constants 
in the denominator; they describe the difference in prop- 
agation velocity for the two different directions. It is 
also quite interesting, that the expression reaches a max- 
imum, if the decay constants are roughly equal. This in- 
dicates a potential resonance in tunnelling circuits, which 
should be detectable by accurate comparisons between 
experiment and theory. Secondly, the denominator con- 
tains the energy differences to the intermediate states in 
the multiple scattering process. This is quite understand- 
able from perturbation theory, where e.g. the interaction 
energies scale with the inverse of the energy difference. 
In this case the intermediate states have to be evaluated 
over the whole energy range. 

Finally, we consider the interaction energy between the 
surface and the tip in the low coupling limit. It has been 
shown recently by an analysis of first order perturbation 
expressions for the tunnelling current and the interaction 
energy, that the two variables should be linear with each 
other. From the first order Greens function we may con- 
struct the density matrix n. The interaction energy is 
then 111: 
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With Gf^^ and Gf^^ 



from Eq. |57| this leads to the 



following result: 
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The energy integration in this case is not trivial, since 
a term containing l//jj, will be singular for 77 — > 0"*" and 
E = Ei or E = Ek- This is due to the missing resistance 
of electrons upon reflection at the surfaces. If a resistance 
term rj is introduced in the same manner as e.g. for the 
contacts F (see Eq. |H)|, then the result converges, and we 
obtain: 
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Physically, we may relate this addition to the fact that 
transiting electrons will lead to charge oscillations in the 
surface layers, depending on the conductance properties 
of the leads. The calculation of the interaction energy 
only involves the computation of the tunnelling matrix 
elements. And, as shown previously, the interaction en- 
ergy will therefore be proportional to the tunnelling cur- 
rent [l3 |. 

To summarize, tunnelling currents and interaction en- 
ergies can be calculated in real space within the non equi- 
librium Greens function formalism based on the separate 
wavefunctions of surface and tip. The zero order expan- 
sion is similar to the traditional Bardeen approach, even 
though we find that each transition carries an individ- 
ual weight. This also implies that resonances between 
surface and tip will play an essential role under specific 
material conditions. The first order expansion describes 
multiple scattering across the tunnelling barrier. In this 
case we can derive a first principles formulation for the 
interaction energy between the two surfaces. 
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